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1.  Summary 

In  this  project,  we  explore  methodologies  to  extract  amplitude  information  from 
empirical  Green  functions  (EGF)  derived  from  ambient  noise  correlations  and  to  map  the 
attenuation  of  the  surface  waves.  Our  objectives  are:  (1)  To  develop  methods  for 
extracting  attenuation  from  ambient  noise;  (2)  To  develop  methods  for  practical 
applications  using  the  Tibetan  Plateau  as  a  test  bed;  (3)  To  develop  preliminary  surface 
wave  attenuation  maps  of  the  Tibetan  Plateau.  Our  approaches  are  to  combine  theoretical 
derivations,  numerical  simulations,  and  practical  considerations.  A  particular  problem  in 
retrieving  amplitudes  from  noise  is  that  seismic  ambient  noise  source  is  not  uniform  and 
it  changes  with  time.  Theoretical  insights  show  that  even  in  the  case  of  incompletely 
diffuse  noise  fields,  we  can  robustly  recover  not  just  travel  times,  but  also  ray  arrival 
amplitudes,  the  ambient  field’s  specific  intensity,  the  strength  and  density  of  its  scatterers 
if  any,  and  most  importantly  attenuation.  We  propose  two  approaches  with  detailed 
formulations:  linear  array  methods  and  more  general  methods  for  2D  station  networks. 
Our  numerical  simulations  validate  that  amplitudes  and  attenuations  can  indeed  be 
extracted  from  noise  correlations  for  a  linear  array  or  for  a  more  general  2D  array.  We 
propose  a  temporal  flattening  procedure  that  is  effective  in  speeding  up  convergence 
while  preserving  relative  amplitudes.  For  real  data,  we  propose  an  “asynchronous” 
temporal  flattening  procedure  that  does  not  require  all  stations  to  have  data  at  the  same 
time.  Tests  on  real  data  suggest  attenuations  are  extracted  that  are  comparable  with  those 
from  earthquakes. 

2.  Introduction 


This  is  the  first  year  of  the  project.  The  objectives  of  this  research  are  to  develop 
methods  for  extracting  attenuation  from  ambient  noise,  to  develop  methods  for  practical 
applications  using  the  Tibetan  Plateau  as  a  test  bed,  and  to  develop  preliminary  surface 
wave  attenuation  maps  of  the  Tibetan  Plateau. 


A  particular  problem  in  retrieving  amplitudes  from  noise  is  that  seismic  ambient  noise 
source  is  not  uniform  and  it  changes  with  time.  Our  approaches  are  to  combine 
theoretical  derivations,  numerical  simulations,  and  practical  considerations.  In  this  report, 
we  show  the  progresses  we’ve  made  on  all  fronts  (theory,  simulation,  and  practice). 

These  progresses  show  not  only  the  great  promise  of  the  methodologies,  but  also  practical 
approaches  for  application  to  real  data. 

3.  Methods  and  procedures 

We  have  described  our  basic  approaches  previously  (Weaver,  2011;  Song  et  al.,  2012). 
Below  is  a  brief  outline  of  the  methodologies  and  procedures.  A  full  examination  of  the 
amplitude  X  of  a  correlation  waveform  between  stations  i  and  j  shows  that  it  takes  the 
form 
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Xff  =  2 1  (m(,  \  xt  -  Xj  |)  exp(-ar  |  x,  -  Xj  |) 


(1) 


where  a.\xj-Xj\  is  the  average  attenuation  between  stations  i  and  j;  co0  is  the  frequency;  c  is 
wavespeed,  and  st  and  sy  are  site  effects  at  the  two  stations. 

B,  is  the  ambient  intensity  in  the  direction  from  i  towards  j  evaluated  at  station  i.  An 
important  and  useful  result  from  this  research  is  that  the  ray  amplitude  X  depends  only  on 
B  in  that  direction  as  an  asymptotically  valid  approximation  (at  the  next  order,  the 
differences  with  respect  to  direction  are  of  order  1%  or  less). 

3.1  The  case  of  a  linear  array 

For  a  linear  array,  along  a  direction  of  several  seismic  stations,  the  amplitudes  X,,  of 

cross  correlation  arrivals  for  j  >  i  is 


(2) 


Formulation  1-1:  path  average.  Eqn  (2)  indicates  that  amplitude  X  from  station  i  (near 
one  end  of  the  array)  to  all  stations  j  in  one  direction  (to  the  other  end  of  the  array) 
depends  only  on  the  intensity  at  station  i  along  the  same  strike  direction.  Thus  the  slope 
of  the  logarithm  of  the  geometrically  corrected  amplitude,  X.sqrt(distance),  as  a  function 
of  distance  would  give  an  estimate  of  the  average  attenuation  (coefficient  a)  along  the 
linear  array. 

Formulation  1-2:  inversion  for  a  linear  array.  The  system  in  eqn  (2)  can  be  linearized 
by  taking  the  logarithm,  as  in  traditional  earthquake-based  attenuation  studies  (e.g.  Yang 
et  al.,  2004).  This  system  is  overdetermined  and  can  be  solved  with  a  reasonable  number 
of  stations.  Taking  N  to  be  the  number  of  stations  along  this  array,  we  have  at  most  N(N- 
l)/2  such  ray-amplitudes  in  one  direction.  Using  one  parameter  per  station  to  describe  the 
local  attenuation,  there  are  up  to  3N  unknowns  (the  site  factor,  the  B  value  along 
direction  n,  and  one  local  attenuation  factor  at  each  station).  Correlations  from  the 
direction  —  n  can  provide  up  to  N(N-l)/2  additional  amplitudes  and  N  additional 
unknowns  B.  Thus  if  amplitudes  along  both  directions  are  available,  one  needs  only  5  to 
7  stations  along  this  line  to  solve  for  all  the  unknowns.  This  is  a  very  general  approach 
without  a  priori  constraints  on  the  noise  intensity. 

Formulation  1-3:  assuming  no  internal  scattering.  The  noise  intensity  cannot  be 
arbitrary.  If  we  assume  the  ambient  noise  sources  are  from  the  far  field  (many  studies 
have  pointed  to  the  oceans  for  the  seismic  bands  we  are  interested  in)  and  further  assume 
that  scattering  is  not  important,  it  can  be  shown  that  the  noise  intensity  B  decays  like 
exp(-2.«. distance),  where  is  a  the  same  medium  attenuation  coefficient  we  are  trying  to 
extract.  Thus  in  this  formulation,  we  need  only  two  Bs  (one  at  each  end  of  the  linear 
array). 


2 

Approved  for  public  release;  distribution  is  unlimited. 


3.2  The  case  of  a  2D  station  network 

Formulation  II-l:  an  expansion  of  the  linear  array  approach.  In  general,  we  have  a 
network  of  stations  whose  distribution  is  irregular  throughout  the  study  region.  However, 
we  may  be  able  to  identify  3  or  more  stations  that  roughly  lie  on  the  same  great  circle 
path.  In  this  case,  we  can  use  amplitude  ratios  to  eliminate  the  unknown  intensity  B  along 
that  direction.  Similar  ideas  are  commonly  used  in  earthquake  studies  to  measure  phase 
velocity  between  two  stations  (the  “two-station”  method)  or  derive  surface  attenuation 
using  amplitude  ratios  between  two  stations  aligned  along  the  great  circle  with  an 
earthquake  (e.g.  Yang  et  al.,  2004). 

Consider  3  stations  (station  1,  2,  3)  on  the  same  great  circle,  the  ratios  of  the 
geometrically-corrected  amplitudes  is  Xi3/Xi2=s3/s2  exp(-a2X2),  X3i/X32=  S1/S2  exp(-aiXi), 
where  s;  is  the  site  factor  for  station  i,  xi  and  X2  are  distance  between  stations  1  and  2  and 
stations  2  and  3,  respectively,  and  oq  and  012  are  average  attenuation  coefficients  between 
stations  1  and  2  and  stations  2  and  3.  If  we  assume  no  internal  scattering,  such  that 
intensity  B  decays  like  exp(-2a.distance)  as  discussed  above,  we  can  use  two  other 
amplitude  ratios:  X23/Xi2=s3/si  exp(-aiXi-a2X2),  and  X2i/X32=si/s3  exp(-aiXi-ot2X2).  The 
logarithm  of  such  an  amplitude  ratio  is  linear  with  distance.  Once  enough  amplitude 
ratios  are  measured  along  paths  of  different  azimuths  and  locations  to  provide  sufficient 
coverage,  it  is  straightforward  to  discretize  the  region  and  construct  a  tomographic 
inversion  for  the  attenuation  coefficients  for  the  entire  region. 

Formulation  II-2:  the  most  general  approach.  The  linear  array  approach  relies  on 
identifying  a  linear  array  of  sufficient  length.  How  might  one  use  sets  of  seismic  stations 
that  are  not  aligned  along  a  linear  array?  How  quickly  ought  one  permit  the  diffuse 
intensity  B  implicit  in  the  above  formulation  to  vary  in  direction  and  from  place  to  place? 
The  diffuse  intensity  B,  in  particular  how  it  is  permitted  to  vary  in  space  and  direction, 
cannot  be  arbitrary.  The  noise  field  ought  to  obey  a  Radiative  Transfer  Equation  (RTE) 
(e.g.,  Turner  and  Weaver  1994),  which  involves  noise  intensity,  noise  sources,  scattering 
strength,  and  attenuation.  In  this  project,  we  will  explore  this  situation  in  theory  and  in 
practice. 

4.  Results 

4.1.  Numerical  simulations:  Tests  on  linear  arrays 

We  reported  previously  that  attenuation  can  be  extracted  along  a  linear  array  using 
simulated  data  with  an  azimuthally  varying  noise  source  (Weaver,  2011;  Song  et  al., 
2012).  The  simulations  included  a  simple  case  of  a  uniformly  spaced  linear  array  of 
receivers,  homogeneous  attenuation  and  wave  speed  and  mild  directionality  to  the  noise 
field.  Numerical  simulations  have  now  been  extended  to  a  broader  class  of  systems,  with 
spatially  varying  attenuation,  highly  directional  noise  intensities,  and  irregular  line  arrays. 
This  is  described  in  detail  in  Weaver  (2013).  However,  the  random  noise  sources  do  not 
change  in  time.  For  these  simulations,  it  is  shown  that  accurate  attenuations  and  site 
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factors  can  be  retrieved.  Furthermore,  errors  in  attenuation  measurements  can  also  be 
determined  (Fig.  1). 


Fig.  1.  Simulation  ‘A’.  Snapshot  of  the  broad-band  wavefield  in  a  271  x  271  array, 
which  seismically  scales  to  800  km,  before  bandpass  filtering.  The  strong  absorption 
on  the  edges  and  anisotropic  annulus  of  sources  are  evident.  The  horizontal  line 
indicates  the  line  along  which  six  stations  were  uniformly  placed  in  the  simulation. 
Other  simulations  had  irregular  station  spacings  and  a  longer  line  array  extending  a  bit 
further  to  the  right.  (Weaver,  2013) 


(a)  (b) 


Lcftgomg  Amplitudes 


Rightgoing  Amplitudes 


Fig.  2.  (a)  The  amplitudes  X  of  the  arrival  wave  packets  in  the  rightgoing 
correlation  waveforms  from  simulation  ‘A’  are  plotted  versus  interstation  distance. 

One-sigma  error  bars  are  taken  from  the  last  column  of  the  table,  SlogX  =  SX/X. 
Amplitudes  corresponding  to  the  same  pseudo-source  fit  well  to  straight  lines.  The 
observed  slopes  correspond  well  to  the  expected  0.00777  nepers  per  mesh  spacing,  (b) 
The  leftgoing  amplitudes  in  simulation  ‘A  ’  have  larger  error  bars,  corresponding  to  the 
relatively  weaker  ambient  noise  intensity  coming  from  the  right. 
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4.2  Numerical  simulations:  “Earth-like”  temporally  variable  noise  field 
1)  A  ID  array  case 

To  examine  our  ability  to  extract  attenuation  from  ambient  noise  methods,  we  impose 
Earth-like  temporally  variable  noise  intensity  on  simulated  diffuse  fields.  The  setup  for 
our  simulation  is  shown  in  Fig.  3.  The  source  intensity  includes  smoothly  varying 
intensity  as  well  as  a  strong  “local”  intensity  (as  indicated  from  the  shade  on  the  source 
“ring”).  The  study  region  is  divided  into  four  parts  inside  the  source  “ring”.  10  receiver 
stations  with  a  spacing  distance  of  85  km  were  arranged  in  a  linear  array,  where  5 
receiver  stations  are  in  the  lower  left  region  and  the  other  5  receiver  stations  are  located 
in  the  upper  right  region.  An  attenuation  model  with  two  different  attenuation 
coefficients,  ax  =  0.00259  in  the  lower  left  region  and  a2  =  0.00388  in  the  upper  right 
region,  was  used.  The  site  responses  of  all  stations  in  the  experiment  are  set  to  1 . 

We  add  the  seismic  data  from  a  station  in  the  USArray  into  our  simulation  data  to 
simulate  the  natural  seismic  noise.  The  seismic  noise  is  from  a  4-year  continuous  seismic 
waveform  after  clipping  earthquakes  and  spikes.  We  then  multiply  this  seismic  noise  with 
the  synthetic  noise  field  point  by  point  to  generate  the  synthetic  noise  field  with 
temporally  varying  strength.  We  call  the  original  simulation  data  “original”  data  and  the 
data  with  seismic  noise  “raw”  data  in  our  discussion  below.  The  sample  rate  of  the  data  is 
3  points  per  second,  and  a  total  of  about  1  year  of  data  are  generated  at  each  station. 

In  order  to  retain  rapid  convergence  amidst  time-varying  noise  intensity,  we  use  a 
“temporal  flattening”  procedure  that  we  proposed  earlier  (Weaver,  2013;  Song  et  ah, 
2012).  Each  station's  band-limited  signal  is  normalized  by  a  running  average  of  the  total 
band-limited  array  energy.  Thus  each  station  is  treated  equally,  and  relative  amplitudes 
are  preserved. 


Fig.  3.  Simulation  ‘B\  (A)  Distribution  of  source 
intensity  and  receiver  locations.  The  Gaussian  noise 
source  is  highly  non-uniform  as  indicated  by  the  shades 
along  the  circumference.  We  have  also  added  temporal 
variation  of  the  noise  source,  which  mimics  temporal 
variation  of  a  seismic  station  on  Earth.  The  attenuation 
inputs  include  four  different  values  at  the  four 
quadrants. 
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Fig.  4.  Amplitudes  extracted  from  Green’s  Functions  with  different  pre-processing 
procedures  for  Simulation  ‘B\  The  lines  are  true  values  from  the  inputs,  (a)  from 
Station  1  to  other  stations,  (b)  from  Station  10  to  other  stations. 

The  results  are  shown  in  Fig.  4.  We  see  that  the  amplitude  decays  of  the  “raw”  data  are 
similar  to  that  of  the  “original”.  But  when  the  signal  to  noise  ratios  of  the  Green’s 
Functions  is  low  (Fig.  4b),  the  amplitudes  are  not  stable.  The  amplitudes  from  the  one-bit 
processing  are  apparently  not  correct;  the  behavior  is  not  linear.  The  amplitude  decays 
from  the  flattening  procedure  are  very  close  to  the  true  values. 

2)  A  2D  array  case 

The  set  up  for  this  simulation  (Simulation  ‘C’)  is  the  same  as  above,  but  the  wavefield  is 
recorded  at  a  2D  array  (10x10  or  a  total  of  100  stations  with  uniform  spacing).  The  input 
attenuation  of  the  medium  is  not  uniform:  it  has  four  different  values  at  the  four 
quadrants  of  the  array.  Following  Formulation  II- 1,  we  identify  linear  arrays  within  an 
azimuth  of  10  degrees  for  each  station.  We  then  form  amplitude  ratios  for  stations  along 
each  linear  array.  We  obtain  a  total  of  about  5600  independent  amplitude  ratios.  These 
amplitude  ratios  are  then  used  to  invert  for  2D  attenuation  maps  and  site  factors.  The 
inversion  with  4  cells  (as  in  the  input  model)  recovers  perfectly  the  input  model  and  the 
inversion  with  36  cells  recovers  quite  well  the  input  model  except  cells  at  the  comers 
which  are  not  constrained  by  the  amplitude  ratios  (Fig.  5). 
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Fig.  5.  Extraction  of  attenuation  from  2D  array  with  simulated  data  (Simulation 

‘C’).  The  input  model  has  four  different  values  at  the  four  quadrants,  respectively.  The 
noise  source  varies  in  strength  azimuthally.  (A)  The  input  Qs  are  fully  recovered  using 
a  4-cell  parameterization.  (B)  The  input  Qs  are  reasonably  recovered  in  a  36-cell 
inversion,  except  at  the  four  corners,  which  are  not  constrained  by  the  method. 


4.3.  Tests  on  real  data 

We  have  conducted  various  tests  on  amplitude  extracted  from  real  data.  Here  we  show 
some  test  results  using  stations  in  the  western  U.S.  (Fig.  6).  We  dealt  with  a  number  of 
practical  issues,  including  energetic  arrivals,  not  all  stations  have  signals,  and  energetic 
arrivals  may  appear  at  different  times  for  different  stations  (Fig.  7),  by  a  new  proposed 
“asynchronous”  temporal  flattening  (ATF)  procedure  (Fig.  8);  this  led  to  an  issue  of 
window  length  for  the  flattening  procedure  (Fig.  9).  The  amplitude  attenuations  extracted 
from  the  ambient  noise  seem  compatible  with  those  from  earthquakes  (Fig.  10). 

Because  not  all  stations  have  signals  and  energetic  arrivals  may  appear  at  different  times 
for  different  stations,  if  we  require  all  stations  have  the  exactly  the  same  overlapping  time 
segments,  we  will  discard  a  lot  of  data,  making  the  convergence  of  the  Green’s  Functions 
difficult.  To  address  this  issue,  we  proposed  a  new  “asynchronous”  temporal  flattening 
(ATF)  procedure,  as  follows. 

1)  Window  out  energetic  signals  for  each  individual  trace.  Keep  track  of  points  that 
are  kept  and  points  that  are  windowed  out. 

2)  Select  a  window  length  for  ATF. 

3)  Calculate  the  RMS  amplitude  of  ALL  traces  for  that  time  window.  Only  points 
that  have  been  kept  after  windowing  are  used. 

4)  Divide  each  trace  within  the  time  window  by  the  RMS. 
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Fig.  6.  Map  of  stations  in  the  western  U.S.  used  in  our  tests.  Indicated  also  are  two 
earthquakes  used  to  compare  with  ambient  noise  results. 
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Fig.  7.  Data  preprocessing.  For  a  given  trace,  energetic  signals  are  first  windowed  out. 
A  flag  trace  is  constructed  to  keep  track  of  points  that  are  kept  (values  of  1)  or 
windowed  out  (values  of  0).  The  windowed  trace  is  then  “flattened”  using  the  average 
root-mean-square  amplitude  over  a  period  of  time  (e.g.  hours)  using  all  stations  and 
data  points  that  have  not  been  windowed  out. 


Fig.  8.  (A)  A  raw  empirical  Greens  function  (EGF)  is  constructed  from  the  cross¬ 
correlation  of  “flattened”  traces.  (B)  Cross-correlation  is  constructed  using  the  flag 
traces,  which  indicates  the  number  of  data  points  that  have  been  used  in  the  EGF 
for  a  given  time  lapse.  (C)  The  raw  EGF  is  normalized  by  the  number  of  data  points 
in  the  flag  CC  to  yield  the  normalized  EGF. 
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Fig.  9.  Tests  on  the  windows  used  in  the  temporal  flattening  procedure.  Shown  are 
results  for  window  length  of  2  hours  to  24  hours.  The  traces  and  amplitudes  are 
identical  and  thus  our  window  length  in  the  flattening  procedure  can  be  quite  wide  (1 
day). 
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Fig.  10.  Comparison  of  amplitude  decays  obtained  from  earthquakes  (red  dots)  and 
ambient  noise  (blue  crosses)  at  two  periods  (10  and  20  s).  Left  two  panels  are  for  the 

western  group  of  stations  and  the  event  in  2007  (Fig.  4)  and  right  panels  are  for  the 
eastern  group  of  stations  and  the  event  in  2009.  Approximately  18  months  of  ambient 
noise  data  are  used  in  calculating  the  EGFs. 

5.  Conclusions 


1)  Based  on  theoretical  considerations,  we  formulated  several  approaches  for 
extraction  of  amplitudes  from  ambient  noise.  The  formulations  take  into  account 
non-uniform  and  anisotropic  source  and  internal  scattering. 

2)  Numerical  simulations  validate  that  amplitudes  and  attenuations  can  indeed  be 
extracted  from  noise  correlations  for  a  linear  array  or  for  a  more  general  2D  array. 

3)  We  propose  a  temporal  flattening  procedure  that  is  effective  in  speeding  up 
convergence  while  preserving  relative  amplitudes.  For  real  data,  we  propose  an 
“asynchronous”  temporal  flattening  procedure  that  does  not  require  all  stations  to 
have  data  at  the  same  time.  The  flattening  procedure  can  have  a  long  window 
length  (1  day). 

4)  Tests  on  real  data  extract  attenuations  that  are  comparable  with  those  from  earthquakes. 
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